function [W, H] = KLNMF(V, r, niter, varargin)

[n,m] = size(V);

W = rand(n,r);
H = rand(r,m);
    
[W H] = Popt(varargin, 'W0', W, 'H0', H);

    
myeps = 1e-6;

% initial reconstruction
R = W*H;

% preallocate matrix of ones
Onm = ones(n,m);

for t = 1:niter

    W = W .* ( ((V./R)*H') ./ max(Onm*H', myeps) );  
    % update reconstruction
    R = W*H;
    
    H = H .* ( (W'*(V./R)) ./ max(W'*Onm, myeps) );
    R = W*H;
    

end


